#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

#include "Vector/vector.03.h"
#include "Color/color.01.h"
#include "Timer/timer.01.h"
#include <vector>
#include <deque> 

using namespace std;


class ICs { 
	public:
	
	typedef unsigned int uint;
	typedef long unsigned int luint;
	typedef double real;
	
	struct P {
		Vector<real> R;
		Vector<real> V;
		int type;
	};
	vector<P> p;
	
	struct PW {
		double e;
		double s; 
		double rc;
	};
	PW pw[10][10];

	struct PT {		
		double m;
		Color color;	
	};
	vector<PT> pt;		
	
	class Grid{
		public:
			int xn, yn, zn;
			double xs, ys, zs;
		
			vector< vector<P*> > node;
			
			inline vector<P*> &operator()(int xi, int yi, int zi) {
				return node[xi + xn*(yi + yn*zi)];
			}
	
	} grid;
	
	
	double kfac;
	
	Timer trf;		
	
	double lU; // length unit
	double mU; // mass unit
	double tU; // time unit
	double eU; // energy unit
	double vU; // velocity unit
	
	double r2_overlap;
	
	int p_size;
	
	ICs():
	lU(1), mU(1), tU(1), 
	eU(mU*(lU*lU)/(tU*tU)), vU(lU/tU),
	kfac(0.977),
	r2_overlap(0.95*0.95)
	{
	}
	
	void addHCPBox(int type1, int type2, Vector<double> R, Vector<int> D, Vector<double> S, Vector<double> V,  double t) {
			
		double dr = pw[type1][type2].s*kfac;
			
		for (int i=0; i<D.x; i++)
		for (int j=0; j<D.y; j++) 
			for (int k=0; k<D.z; k++){
					
				P pi;	
					
				pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
				pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
				pi.R.z = sqrt(6.0)/3.0*dr*k;
					
				//pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));
				pi.V *=0;
				
				//pi.R.x-= 0.5*S.x*dr;
				//pi.R.y-= 0.5*S.y*0.5*sqrt(3)*dr;
				//pi.R.z-= 0.5*S.z*sqrt(6.0)/3.0*dr;
				
				pi.R.x*= S.x;
				pi.R.y*= S.y;
				pi.R.z*= S.z;
				
				
				pi.R  += R;
				pi.V  += V;
				
				pi.type =  (i==0 || j==0 || k==0 || i==D.x-1 || j==D.y-1 || k==D.z-1) ? type2 : type1;

				
				
				p.push_back(pi);	
					
			}
		p_size = p.size();
	}
	
	
	inline void find_overlaps(vector<P*> &pi) {
	
	
		for (int i=0;i<pi.size()-1;i++) if (pi[i]->type!=-1)
		for (int j=i+1;j<pi.size();j++) {
						
			real dx = pi[i]->R.x - pi[j]->R.x; // 1
			real dy = pi[i]->R.y - pi[j]->R.y; // 1
			real dz = pi[i]->R.z - pi[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 < r2_overlap) {
				if (pi[i]->type==0) pi[j]->type=0;
				pi[i]->type =-1;
				continue;
			}
		}
	}


	inline void find_overlaps(vector<P*> &pi,  vector<P*> &pj) {
					
		for (int i=0;i<pi.size();i++) if (pi[i]->type!=-1)
		for (int j=0;j<pj.size();j++) {
						
			real dx = pi[i]->R.x - pj[j]->R.x; // 1
			real dy = pi[i]->R.y - pj[j]->R.y; // 1
			real dz = pi[i]->R.z - pj[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 < r2_overlap){
				if (pi[i]->type==0) pj[j]->type=0;
				pi[i]->type =-1;
				continue;
			}
				
			
		}
					
	}
	
	void remove_overlaps(){

		trf.set();
		
		grid.node.resize(grid.xn*grid.yn*grid.zn);
	
		for (int i=0; i<grid.node.size(); i++) grid.node[i].resize(0);
		
		
		double grid_05lenx = 0.5*grid.xn*grid.xs;
		double grid_05leny = 0.5*grid.yn*grid.ys;
		int k=0;
		p_size = p.size();
		for(int i=0; i<p_size; i++){
 			if (fabs(p[i].R.x) > grid_05lenx) { printf( "Warning grid_05lenx = %g < |%g = p[%i].R.x| we better stop!\n",i, grid_05lenx, p[i].R.x); exit(0); }
 			if (fabs(p[i].R.y) > grid_05leny) { printf( "Warning grid_05leny = %g < |%g = p[%i].R.y| we better stop!\n",i, grid_05leny, p[i].R.y); exit(0); }
 			
			grid(int((p[i].R.x+grid_05lenx)/grid.xs), int((p[i].R.y+grid_05leny)/grid.ys), int(p[i].R.z/grid.zs)).push_back(&p[i]);
			
		}
		
		//printf("remove_overlaps(): grid has been saccessfully assigned\n");
		
		
		
		for (int xi=0; xi<grid.xn; xi++) 
		for (int yi=0; yi<grid.yn; yi++)
		for (int zi=0; zi<grid.zn; zi++) {
			
			vector<P*> &pi = grid(xi, yi, zi);
			
			if (pi.size()==0) continue;
							
			find_overlaps(pi);
			
			if (0 < xi && 0 < yi && zi+1<grid.zn)					find_overlaps(pi, grid(xi-1,yi-1,zi+1)); // 12
			if (0 < xi && zi+1<grid.zn)							   find_overlaps(pi, grid(xi-1,yi  ,zi+1)); // 1			
			if (0 < xi && yi+1<grid.yn)								find_overlaps(pi, grid(xi-1,yi+1,zi  )); // 13
			if (0 < xi && yi+1<grid.yn && zi+1<grid.zn)			find_overlaps(pi, grid(xi-1,yi+1,zi+1)); // 6
			
			if (0 < yi && zi+1<grid.zn)								find_overlaps(pi, grid(xi  ,yi-1,zi+1)); // 4
			if (zi+1<grid.zn)												find_overlaps(pi, grid(xi  ,yi  ,zi+1)); // 2
			if (yi+1<grid.yn && zi+1<grid.zn)						find_overlaps(pi, grid(xi  ,yi+1,zi+1)); // 5
			if (yi+1<grid.yn)												find_overlaps(pi, grid(xi  ,yi+1,zi  )); // 8
			
			if (xi+1<grid.xn && 0 < yi && zi+1<grid.zn)			find_overlaps(pi, grid(xi+1,yi-1,zi+1)); // 11
			if (xi+1<grid.xn && zi+1<grid.zn)						find_overlaps(pi, grid(xi+1,yi  ,zi+1)); // 3
			if (xi+1<grid.xn )											find_overlaps(pi, grid(xi+1,yi  ,zi  )); // 10
			if (xi+1<grid.xn && yi+1<grid.yn && zi+1<grid.zn)	find_overlaps(pi, grid(xi+1,yi+1,zi+1)); // 7
			if (xi+1<grid.xn && yi+1<grid.yn)						find_overlaps(pi, grid(xi+1,yi+1,zi  )); // 9
							
		}
	
		trf.get();
	
		for (int i=0; i<p_size; i++) 
			if (p[i].type == -1) p[i--] = p[--p_size];	
		
	}
	
	void save(string filename) {
			
		FILE *fp = fopen (filename.c_str() , "wb");
		
		struct MEM { Vector<float> R; int type;} *mem;
	
		mem = new MEM[p_size];
		for (int i=0; i<p_size; i++) {
			
			if (p[i].type!=0 && p[i].type!=1) printf("ICs::save() Warining!\n\tp[%i].type!=0 && p[%i].type!=1  = %i\n", i,i, p[i].type);
			
			mem[i].R = p[i].R;
			mem[i].type = p[i].type;
		}
		
		fwrite(mem, sizeof(MEM), p_size, fp);
		fclose(fp);
		delete [] mem;
	
	}
	
	void print(){
		
		printf("ICs::print()\n");
		
		printf("\tremoved particles: %i\n", p.size()-p_size);
		printf("\tfinding of overlaps: %g sec \n", trf.sin);
		printf("\tnumber of particle: %i\n", p_size);
		//printf("finding of overlaps: %g sec \n", trf.sin);
		//double g_memsize = ;
		printf("\tgrid.node[0]: %i byte\n", sizeof(grid.node[0]));
		printf("\tsizeof(grid.node[0]))*grid.node.size() = %g Mbyte\n", sizeof(grid.node[0])*(float)grid.node.size()/1024/1024);
		printf("\tgrid size: %g MByte\n", (sizeof(grid.node[0])*(float)grid.node.size() +  sizeof(P*)*p.size())/1024/1024);
		printf("\n");
		//double p_memsize = sizeof(p[0])*double(p.size());
		
	}
	
	
	void make_WTC() {
	
		grid.xn = 2*64;
		grid.yn = 2*64;
		grid.zn = 300;
	
		grid.xs = 3;
		grid.ys = 3;
		grid.zs = 3;
	
		pt.resize(2);
	
		pt[0].m = 1;
		pt[1].m = 1;
	
		pw[0][0].s  = 1;
		pw[0][0].rc = 3;//SQR(2);
	
		pw[1][1].s  = 1;
		pw[1][1].rc = 3;//SQR(2);
	
		pw[0][1].s  = pw[1][0].s = 1;
		pw[0][1].rc = pw[1][0].rc = 3;//SQR(2);	

		lU = 1.0/6; // meters
		tU = 1; //sec
		vU = lU/tU;
	
		
		// core columns
		addHCPBox(0, 1, dVector( ( int(-20/lU)-0  )*kfac, (int(+12/lU/(.5*sqrt(3)))  )*(.5*sqrt(3))*kfac, 0), Vector<int>(12,8,350), dVector(1,1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( ( int(-20/lU)-0  )*kfac, (int(-12/lU/(.5*sqrt(3)))-8)*(.5*sqrt(3))*kfac, 0), Vector<int>(12,8,350), dVector(1,1,1), dVector(0,0,0), 0);
		
		addHCPBox(0, 1, dVector( ( int(+20/lU)-12 )*kfac, (int(+12/lU/(.5*sqrt(3)))  )*(.5*sqrt(3))*kfac, 0), Vector<int>(12,8,350), dVector(1,1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( ( int(+20/lU)-12 )*kfac, (int(-12/lU/(.5*sqrt(3)))-8)*(.5*sqrt(3))*kfac, 0), Vector<int>(12,8,350), dVector(1,1,1), dVector(0,0,0), 0);
		
		
		// facade columns
		addHCPBox(0, 1, dVector( ( int(-31/lU)-0 )*kfac, (int(+31/lU/(.5*sqrt(3)))-7)*(.5*sqrt(3))*kfac, 0), Vector<int>(8,8,350), dVector(1,1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( ( int(-31/lU)-0 )*kfac, (int(-31/lU/(.5*sqrt(3)))-1)*(.5*sqrt(3))*kfac, 0), Vector<int>(8,8,350), dVector(1,1,1), dVector(0,0,0), 0);
		
		addHCPBox(0, 1, dVector( ( int(+31/lU)-8 )*kfac, (int(+31/lU/(.5*sqrt(3)))-7)*(.5*sqrt(3))*kfac, 0), Vector<int>(8,8,350), dVector(1,1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( ( int(+31/lU)-8 )*kfac, (int(-31/lU/(.5*sqrt(3)))-1)*(.5*sqrt(3))*kfac, 0), Vector<int>(8,8,350), dVector(1,1,1), dVector(0,0,0), 0);

		for (int fh=36; fh<350; fh+=36) {
		
		// internal beems
		addHCPBox(0, 1, dVector( int(-31/lU)*kfac,  (int(+12/lU/(.5*sqrt(3)))   + 6)*(.5*sqrt(3))*kfac, fh*sqrt(6.0)/3.0*kfac), Vector<int>(int(+62/lU),8,8), dVector(1,1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( int(-31/lU)*kfac,  (int(-12/lU/(.5*sqrt(3)))-8 - 6)*(.5*sqrt(3))*kfac, fh*sqrt(6.0)/3.0*kfac), Vector<int>(int(+62/lU),8,8), dVector(1,1,1), dVector(0,0,0), 0);
		
		addHCPBox(0, 1, dVector( ( int(+20/lU) -2  ) *kfac, ( int(-31/lU/(.5*sqrt(3)))-1 )*(.5*sqrt(3))*kfac,  fh*sqrt(6.0)/3.0*kfac),  Vector<int>(8,int(+62/lU/(.5*sqrt(3))),8), dVector(1,+1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( ( int(-20/lU) +2-8) *kfac, ( int(-31/lU/(.5*sqrt(3)))-1 )*(.5*sqrt(3))*kfac,  fh*sqrt(6.0)/3.0*kfac),  Vector<int>(8,int(+62/lU/(.5*sqrt(3))),8), dVector(1,+1,1), dVector(0,0,0), 0);		
		
		// facade beems
		addHCPBox(0, 1, dVector( int(-31/lU)*kfac,  (int(+31/lU/(.5*sqrt(3)))-7)*(.5*sqrt(3))*kfac, fh*sqrt(6.0)/3.0*kfac), Vector<int>(int(+62/lU),8,8), dVector(1,1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( int(-31/lU)*kfac,  (int(-31/lU/(.5*sqrt(3)))-1)*(.5*sqrt(3))*kfac, fh*sqrt(6.0)/3.0*kfac), Vector<int>(int(+62/lU),8,8), dVector(1,1,1), dVector(0,0,0), 0);
		
		addHCPBox(0, 1, dVector( ( int(+31/lU) -8) *kfac, ( int(-31/lU/(.5*sqrt(3)))-1 )*(.5*sqrt(3))*kfac,  fh*sqrt(6.0)/3.0*kfac),  Vector<int>(8,int(+62/lU/(.5*sqrt(3))),8), dVector(1,+1,1), dVector(0,0,0), 0);
		addHCPBox(0, 1, dVector( ( int(-31/lU) -0) *kfac, ( int(-31/lU/(.5*sqrt(3)))-1 )*(.5*sqrt(3))*kfac,  fh*sqrt(6.0)/3.0*kfac),  Vector<int>(8,int(+62/lU/(.5*sqrt(3))),8), dVector(1,+1,1), dVector(0,0,0), 0);		

		}
		
		
		
		
		remove_overlaps();
		print();
		save("ics.dat");	
		
		grid.node.resize(0);
		//exit(0);
	
		mU = 1;//3.87e6*110/p.size(); //kg
		eU = mU*SQR(lU)/SQR(tU); // j

		
		pw[0][0].e  = SQR(tU)/SQR(lU); 
		pw[1][1].e  = SQR(tU)/SQR(lU);
		pw[0][1].e  = pw[1][0].e = SQR(tU)/SQR(lU);

	}
	
};
